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The present study is aimed at development of a dynamic model for the beta-type Stirling engines with 
rhombic-drive mechanism. Dynamic simulation of the engine is carried out under different operating 
and geometrical conditions by connecting the dynamic model to the existing thermodynamic model 
[23]. In the connection, the instantaneous gas forces exerted on the piston and the displacer are 
determined from the gas pressures in the expansion and the compression chambers with the help of the 
thermodynamic model. Once the gas forces are obtained, the dynamic model is used to determine the 
positions, velocities, and accelerations of the components of the engine at the consecutive time step, and 
then the gas pressures in the expansion and the compression chambers can be updated by the ther¬ 
modynamic model. In this manner, transient variation in rotational speed of the engine during the hot- 
start period and the performance curves of the engine indicating the dependence of the power output 
and the thermal efficiency on the rotation speed can be predicted. In the present study, the effects of the 
major geometrical and operating parameters on the steady-state rotational speed, maximum power 
output and thermal efficiency have been evaluated thoroughly in a comprehensive parametric study. 

© 2011 Elsevier Ltd. All rights reserved. 


1. Introduction 

Concentrating solar power systems with Stirling engines (CSP- 
Stirling) have been treated as one of the most promising renewable 
energy technologies and received great attention from the related 
researchers in the past several years [1—3]. The capacity of a typical 
CSP-Stirling units ranges from 3 to 25 kW. In general, the CSP- 
Stirling system is equipped with a solar energy absorber to 
receive the solar radiation and a Stirling engine to convert the 
received solar energy to the electrical power output. Theoretically, 
the performance of the CSP-Stirling system for generating electrical 
power is strongly dependent on the energy conversion efficiencies 
of the solar energy absorber and the Stirling engine. According to 
the earlier review of Ross [4], among all different sizes, the 
domestic-scale Stirling electric power generator [5] is particularly 
of great market potential. 

Stirling engines operate over a closed, regenerative thermody¬ 
namic gas cycle, with cyclic compression and expansion of the 
working gas in different chambers at different temperature levels. 
The working gas in the Stirling engines might be air, nitrogen, 
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helium, or hydrogen. In terms of the arrangement of piston and 
displacer, the Stirling engines are classified into three types of 
configuration, namely alpha-, beta- and gamma-type [6]. In the 
past several decades, a number of experimental studies of the 
performance of different types of Stirling engines have been con¬ 
ducted [7—11]. A detailed literature review regarding the Stirling 
engines design was provided by Kongtragool and Wongwises [12]. 

On the other hand, numerical simulation based on theoretical 
model is also regarded as a useful tool toward understanding of the 
physical transport phenomena in the Stirling engines that might not 
be easily obtained by experiments. Schmidt [13] firstly proposed an 
isothermal model which has been used widely due to its mathe¬ 
matical simplicity. The main assumption of Schimdt’s theory is that 
the chambers inside the engine are of the same pressure and each 
chamber is fixed at a constant but distinct temperature. Lately, 
Finkelstein [14] performed a non-isothermal analysis and evaluated 
the heat transfer coefficient and the temperature variation in the 
working space. Makhkamov and Ingham [15] divided the internal 
circuit of the Stirling engine into five chambers, and solved the 
governing ordinary differential equations for the energy and mass 
conservation. Mechanical losses in the construction of the engine 
have been determined and the computational results show that the 
mechanical losses for this particular design of the Stirling engine 
may be up to 50% of the indicated power of the engine. Recently, the 
















162 


C.-H. Cheng, Y.-J. Yu /Renewable Energy 37 (2012) 161—173 


Nomenclature 

R d 

offset distance from gear centers to joints between 




rhomboid and gears (m) 

a 

acceleration (m/s 2 ) 

Rn, R t2 thermal resistances of expansion and compression 

Qoad 

load coefficient 


chambers (°C/kW) 

dt 

time step (s) 

t 

time (s) 

F, 

friction force at contact surface between piston and 

tp 

cycle period (s) 


cylinder (N) 

T h 

heat source temperature (K) 

f 2 

friction force at contact surface between displace rod 

Tc 

heat sink temperature (K) 


and piston (N) 

w 

steady-state output power (W) 

G 

regenerative channel gap (m) 

X,Y 

internal forces in x- and y-direction (N) 

h 

flywheel moment of inertia (kg m/s 2 ) 


rectangular coordinates 

k 

length of displacer (m) 

y 

displacement (m) 

M 

moment (N m) 

Z\ to z 4 geometrical constants, Eqs. (2a)—(2d) 

m\ 

mass of displacer assembly (displacer, displacer rod, 




and link 7-8) (kg) 

Greek Symbols 

m 2 

mass of piston (kg) 

V 

thermal efficiency 

m 3 

mass of left gear (kg) 

0 

angle (rad) 

m 4 

mass of right gear (kg) 

lx) 

instantaneous rotational speed (rad/s) 

m 6 

mass of link 6-10 and piston connecting rod (kg) 

b)ave 

cycle-average rotational speed (rad/s) 

m 7 

mass of link 5-6 (kg) 

a 

angular acceleration (rad/s 2 ) 

m 8 

mass of link 5-7 (kg) 

0 

phase angle (degree) 

m 9 

mass of link 8-9 (kg) 

T 

torque (N m) 

mio 

mass of link 9-10 (kg) 

V 

velocity (m/s) 

™air 

mass of working air (kg) 



P 

pressure (kPa) 

Subscripts 

Qin 

input heat transfer rate (W) 

initial 

initial value 

Q 

crank angle (rad) 

d 

displacer 

r 

radius (m) 

load 

external load 



P 

piston 


progress in the development of numerical models for the Stirling 
engines has been facilitated significantly [16-21]. 

Rhombic-drive mechanism was proposed by the Philips 
company [22], and has been commonly adopted in the |3-type 
Stirling engines. In the rhombic-drive mechanism, one rigid rod is 
installed connecting the piston to the top corner of a jointed 
rhomboid, and another rigid rod connecting the displacer to the 
bottom corner of the rhomboid. In addition, two symmetric gears of 
equal diameter are connected to the right and the left corners of the 
rhomboid fixed on the gears at an offset distance from gears center. 
The two gears are in contact and rotate in opposite directions, one 
of which is connected to the flywheel. When gas pressure is applied 
to the piston, the top corner of the rhomboid is pushed downward; 
the rhomboid is flattened in the direction of the piston axis; and 
then it pushes on the gear wheels and causes them to rotate. As the 
gear wheels rotate, the rhomboid experiences its change of shape 
by driving its top and bottom corners to move upwards and 
downwards, respectively. In this manner, the piston and the dis¬ 
placer are moved back-and-forth in the cylinder in a perfectly 
coaxial way. The major advantages of the rhombic-drive mecha¬ 
nism are that the coaxial movement of the piston and the displacer 
is quiet and requires no serious lubrication. Recently, Cheng and Yu 
[23] presented a thermodynamic model for the beta-type Stirling 
engine with the rhombic-drive mechanism. Periodic variation of 
pressures, volumes, temperatures, masses, and heat transfers in the 
expansion and the compression chambers are predicted by taking 
into account the effectiveness of the regenerative channel as well as 
the thermal resistances of the heating and cooling heads. 

As stated earlier, a number of the thermodynamic models 
[13-21,23] are available for analysis of the performance of the 
Stirling engines; however, the existing models are limited to the 
thermodynamic analysis or mechanical loss in the steady-state 
regime and not valid for the dynamic simulation of the engines in 


the start-up period. Cheng and Yu [24] firstly performed an efficient 
dynamic simulation model for the Stirling engines and investigated 
the transient variation in rotational speed of the engine during the 
hot-start period. Unfortunately, the model is only valid for the cam- 
drive mechanism and cannot be applied for the engines with 
rhombic-drive mechanism; therefore, a dynamic model associated 
with the rhombic-drive engines is still lacking. For an internal 
combustion engine, Giakoumis and Rakopoulos [25] provided 
a numerical model accounting for the transient conditions with the 
slider-crank-driving mechanism. The crankshaft angular 
momentum equilibrium was formulated by taking into account 
instantaneous gas, inertia and friction loads. However, the model 
for the internal combustion engines developed cannot be directly 
applied for the beta-type Stirling engines. 

Under these circumstances, the present study is aimed at devel¬ 
opment of a dynamic model for the beta-type Stirling engines with 
the rhombic-drive mechanism. The dynamic model is built based on 
force and moment balances associated with all the links and 
components of the rhombic-drive mechanism. Dynamic simulation 
of the engine is carried out under different operating and geometrical 
conditions by combining the present dynamic model with the ther¬ 
modynamic model. Effects of the major geometrical and operating 
parameters on the rotational speed and the performance curves of 
the engine are then investigated in a parametric study. 

2. Dynamic model 

2.1. Kinematic analysis 

Fig. 1 shows a typical |3-type Stirling engine with the rhombic- 
drive mechanism and all its parts. The 10 W engine (Model 10ST1) 
is designed and made in the Power Engines and Clean Energy 
Laboratory (PEACE Lab.), National Cheng Rung University. The 
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Fig. 1 . A (3-type Stirling engine with rhombic-drive mechanism and all its parts (Model 10ST1). 


schematic of the engine is illustrated in Fig. 2, where 0 represents the 
crank angle; 85 the angle of the right gear; Rd the offset distance from 
gear centers to the joints between the rhomboid and the gears; l\, I 2 , 
I 2 , and I 4 the lengths of the links; Ldt the length from the linkage l\ to 
the top surface of the displacer; L pt the length from the link l 4 to the 
top surface of the piston; L the distance between the center of the 
two gears; G the size of the regenerative channel; and T 2 the 
radiuses of the displacer and the cylinder, respectively. Basically, the 
schematic is referred to as a planar structure [26], and the piston, 
displacer, and the links of the mechanism are labeled in Fig. 3. In 
practices, the displacer of the engine moves ahead of the piston at 
a fixed phase angle in order to drive the working fluid to the 
expansion chamber to be heated and to the compression chamber to 
be cooled so that the working air experiences expansion and 
compression periodically so as to deliver power through the piston. 

Let q = 8 , the crank angle, and then the angles 0 5 ,0 7 , 6 g, 89 and 6 \o 
are expressed as: 


6 = q 

(la) 

8 5 = Tt-q 

(lb) 

0 7 = cos -1 {Z\ +z 2 cos (q)) 

(lc) 

0 8 = sin" 1 (z 3 +Z 4 COS (q)) + | 

(Id) 

d 9 = n- (sin -1 (z 3 +z 4 cos (<j)) +£) 

(le) 

010 = JT- (cos" 1 (Zj +z 2 cos (q))) 
where z 1 , Z 2 , Z 3 , and Z 4 are defined as 

(If) 
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Fig. 2. Schematic of the beta-type Stirling engine with rhombic-drive mechanism. 


Z\ = L - l\/ 2 l 2 

(2a) 

OJ4 = —CO 


z 2 - Rd/h 

Z3 = L — I4/2I3 

(2b) 

(2c) 

w 5-6 = 

V 

z 2 sin (q)co 

h - (z-\ + Z 2 c os (q)) 2 

z 4 = -Rd/l 3 

(2d) 

W 5 _7 = 

-z 4 sin (q)co 

The rotating motions of the components and the links are rep¬ 
resented by their rotation speeds, which can be determined by 
differentiating Eqs. (la)—(If). That is, 

V 

W9-8 = 

O - (z 3 +Z4COS ( q )) 2 

z 4 sin (q)w 

/ 7 

6J3 = (O 

(3a) 

V 

- (z 3 +Z4COS ( q)) z 
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Fig. 3. Free body diagrams of the components in the rhombic-drive mechanism. 


w 9-10 


-z 2 sin (q)(o 
\/l - (Zi +z 2 cos (q )) 2 



where w is the rotational speed of the engine. If l 2 = I 3 , then #g = #10 
and #1 = 0g, and one has wg = wio and 017 = wg. Next, the rotational 
accelerations are derived by further differentiating the expressions 
of the rotational speeds. The results are as follows: 


= a 

(4a) 

a 4 = -a 

(4b) 


z 2 cos (q)co 2 + z 2 sin (q)a 

“s-e = — 7 - - 

V 1 - (Zi +z 2 cos ( q)) z 

(z 2 sin (q)w) 2 (z-, +z 2 cos (q)) 

T 71 1 ' 5 

l-(Zi+z 2 cos (q)Y 


Z 4 COS (g)w 2 +z 4 sin (g)a 

^ 5-7 =- 7 .= ~ 

yl - (z 3 +Z 4 COS (q )) 2 

(z 4 sin (q)w) 2 (z 3 +Z 4 COS (q)) 
~ l r 7l 1.5 

1 - (z 3 +Z 4 C 0 S (q)) z 


(4c) 


(4d) 
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^9-8 


z 4 cos (g)w 2 +z 4 sin (q)a 


yl-fe +Z 4 COS (q)) : 


(z 4 sin (q)w) 2 (z 3 +z 4 cos (q)) 


1 - (z 3 +z 4 cos (g))‘ 


1.5 


(4e) 


^9-10 


z 2 sin (q)a + z 2 cos (g)w 2 ’ 


+ 


\/l - (Zi +z 2 cos (q)) 2 
(z 2 sin (q)w) 2 (Z! +z 2 cos (q)) 


1 - (zi + z 2 cos (q))' 


1.5 


(4f) 


It is noted that the displacer, the piston, link 6-10, and link 7-8, 
experience no rotation but translation. 

Displacements of the centers of mass of the displacer and the 
piston [y x (q) and y 2 (q)] and the links [y 3 (q). y 4 (q). y 6 -io (q) y 7 -s 
(q), y 5 _ 6 (q), y 5 _ 7 (q), y 9 - 8 (q), and y 9 _ 10 (q)] are derived as: 


The expressions for velocities of the components and the links 
in the y-direction are obtained by differentiating the displacement 
equations, Eqs. (5a)-(5j), with respect to time. The results are 


ki = KdCos (q)q + / 3 sin (d 8 - 

(6a) 

v 2 = R d cos (q)q + / 2 sin (0 7 )0 7 

(6b) 

v 3 = 0 

(6c) 

v 4 = 0 

(6d) 

d6-io = Kdcos (q)q + l 2 cos (d 7 )d 7 

(6e) 

9 7 -8 = RdCos (q)q + / 3 sin (d 8 - |)0 8 

(6f) 


1 

R d cos (q)q + -/ 2 cos (0 7 )d 7 

( 6 g) 

Rd cos (q)q + 7/3 sin 0 8 

( 6 h) 

R d cos (q)q + l/ 3 sin (fl 8 -|) #8 

( 6 i) 

1 

K d cos (q)q + ^cos ( 0 7 ) 0 7 

( 6 j) 


The accelerations in the y-direction are obtained by differenti¬ 
ating the velocity equations. That is 

ax = R d cos (q)q - i? d sin (q)q + 1 3 cos (d 8 ~ 7 ) ®8 


7T 


+ / 3 sin (d 8 -~)d 


8 


yi (<J) = L dt + R d sin (q) - / 3 cos ((? 8 - 7 ) 

(5a) 

a 2 

y 2 (q) = Ut + K d sin (q) + / 2 cos ( 0 7 ) 

(5b) 

a 3 

y 3 (q) = 0 

(5c) 

a 4 

y 4 (q) = 0 

(5d) 

a 6 

ye-io(<?) = «dSin (q) + / 2 sin ( 0 7 ) 

(5e) 

a? 

y 7 _ 8 (q) = i? d sin (q) - / 3 cos ( 0 8 - 

(5f) 


y 5 -e(q) = KdSin(q)+l/ 2 sin( 0 7 ) 

(5g) 

as 

y 5 - 7 (q) = ^dSin (q) - I/ 3 COS (0 8 -|) 

y9-s(<J) = ^dSin (q) - 7/3 cos (d 8 - 1) 

(5h) 

«5 

(5i) 


yg-io(q) = ^sin (q) +l/ 2 sin ( 0 7 ) 

(5j) 



(7a) 

2 -2 

a 2 = R d cos(q)q-R d sm(q)q + l 2 sin(d 7 )0 7 + l 2 cos(d 7 )0 7 (7b) 

(7c) 
(7d) 

a 6 _ 10 = R d cos(q)q-R d sm(q)q z -l 2 sm(d 7 )d 7 + l 2 cos(d 7 )d 7 (7e) 

7T\ k2 


a 7 _ 8 = ^cos (q)q - K d sin (q)q 2 + / 3 cos (tf 8 - 


7P 


+ / 3 sin 0 8 -- (9 


8 


2 1 

1 + 2 


/ 2 sin ( 0 7 )O : 


4 - / 2 cos ( 0 7 )'@ 7 


2 1 

K d cos (q)q - K d sin (q)q + - 


7T' 


/ 3 cos (6>s -7)^ 


+ / 3 sin (0 8 — 2 J 


2 1 


ag-s = fld cos (<?)q - ^d sin (<?)<? + 


/ 3 COS (0 8 




7T 


+ / 3 sin ^ 0 8 - 2)^8 


« 9 -io = ^cos (q)q + l 


- l 2 sin (# 7 )# 7 + / 2 cos (# 7 )# 7 


(7f) 


(7g) 


(7h) 


(7i) 


(7j) 


2.2. Dynamic analysis 

In order to predict the dynamic behavior of the mechanism, the 
free body diagrams for all the components are plotted and shown in 
Fig. 3. The dynamic model is built based on force and moment 
balances associated with all the links and components of the 
rhombic-drive mechanism. It is important to notice that by using 
the same concept of the friction cycle described in Ref. [25], the 
friction moments in the hinge joints, piston sliding surface, and the 
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contact surface of the two gears are included in the model. In 
addition, the effects of the weights of all components are taken into 
account in the model. The force and the moment balance equations 
associated with each free body are written as follows: 

1. Displacer, displacer rod and link 7-8 [Fig. 3(a)]: 

Force balance in x-direction: 

X 8 -X 7 = m 1 a xl (8a) 

Force balance in y-direction: 

-Pe + Y 7 + Y 8 +F 2 -m^g = m^a y1 (8b) 

Moment balance equation: 

(Y s - Y 7 ) l j + Ms - M 7 = 0 (8c) 

2. Piston [Fig. 3(b)]: 

Force balance in x-direction: 

x 2 = -Np + m 2 a x2 (9a) 

Force balance in y-direction: 

-Pc -Y 2 - n'Np - m 2 g = m 2 a y2 (9b) 

Moment balance equation: 

Mp = - M 2 (9c) 

where -P e + Y 7 + Y 8 + F 2 - rrqg = m\a y \ is the assumed friction 
force exerted at the contact surface between the piston and the 
cylinder. 

3. Link 6-10 and piston connecting rod [Fig. 3(c)]: 

Force balance in x-direction: 

X 2 +^10 _ ^6-10 = m 6 a x6-W (10a) 

Force balance in y-direction: 

Y 2 — Y 6 — Y 10 — F 2 — m 6 g = m 6 %6-io (10b) 

Moment balance equation: 

-X 2 (l pt - 0.5/ p ) + 0.5/! (Y 6 - Yio) + M 10 - M 6 - M 2 = 0 (10c) 


4. Link 5-7 [Fig. 3(d)]: 


Force balance in x-direction: 


X 57 +X 7 = m 8 a x5 _ 7 

(1 la) 

Force balance in y-direction: 


Y 57 — Y 7 — m 8 g = m 8 a y5 _ 7 

(lib) 


Moment balance equation: 


/ 3 sin (7T-# 8 )X 7 - Y 7 l 3 cos(7r-0 8 ) 


5. Link 5-6 [Fig. 3(e)]: 

Force balance in x-direction: 

^56-^6 - m 7 a x5-6 (12a) 

Force balance in y-direction: 

y 56 - Y 6 - m 7 g = m 7 a y5 _ 6 (12b) 

Moment balance equation: 

-/ 2 cos ( 0 7 )Xg + Yg/ 2 sin ( 0 7 ) - 2 ^ + +M 56 

= k- 6 a 5-6 (12c) 


6. Link 8-9 [Fig. 3(f)]: 

Force balance in x-direction: 

“(^89+^8) = m 9 a x8-9 (l^a) 

Force balance in y-direction: 

Y 89 - Y 8 - m 9 g = m 9 a y8 _ 9 (13b) 

Moment balance equation: 

-Mn (* 9 )X S + V s l 3 cos <«,) + + M S - M a , 

= h-9 a 8-9 (13c) 


7. Link 9-10 [Fig. 3(g)]: 

Force balance in x-direction: 


910—X 10 = m w d x9 -10 

(14a) 

Force balance in y-direction: 


(Y 9 io - Yio) - rriwg = m io%9-io 

(14b) 

Moment balance equation: 



/ 2 sin (7t- ^io)Xio - Yio/ 2 c°s (7r— /9io) 


8. Right gear [Fig. 3(h)]: 

Force balance in x-direction: 

X 9 io+X 8 9 +X 4 = 0 (15a) 

Force balance in y-direction: 

—(Y 910 - Y 89 ) + Y 34 + Y 4 - m 4 g = 0 (15b) 

Moment balance equation: 

(1^910 + Y 89 )R d cos ( 7 r - 0 5 ) - (X 910 +X 89 ) R^sin (^ - ^ 5 ) 

— Y 34 - + M 910 + M 89 = / 4 a 4 


2 


(15c) 
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9. Left gear [Fig. 3(i)]: 

Force balance in x-direction: 

X 5 6 -X57 -X 3 = 0 ( 16 a) 

Force balance in y-direction: 

^ 56-^57 + ^ 3-0135 = 0 ( 16 b) 

Moment balance equation: 

-(*56 -*57) fidSin ( 6 ) + (Y 56 - Y 57 ) R d c os ( 6 ) 

- Y34- + M 3 - M 57 - M55 - i load - l 3 a 3 ( 16 c) 

where M 3 (= J n a) is the torque induced by the flywheel on the left 
gear, and a 3 is exactly equal to the engine’s flywheel rotational 
acceleration, a, because the left gear is fixed to the flywheel. The 
magnitude of a 3 can be expressed by 

a 3 = ^3 +1 - co^j jdt ( 17 ) 

where n indicates the time step. Combining Eqs. (16c) and 
(17) leads to 

w 3 +1 — w 3 + |^"[ T power — ^load] (18) 

where 

"tpower = — (^56 —^57) fid s i n ($) + (^56 - ^ 57 ) Rd cos ($) 

- Y 34 ^ - M 57 - M 56 ( 19 a) 

^load = Qoad w (19b) 

Note that in the above equation, T p0W er is a combination of the 
effects of the gas power torque, mass inertia torque, and the friction 
torques (M 57 and M 56 ) exerted on the left gear as well. On the other 
hand, Ti oa d denotes the resisting torque resulting from the external 
load. The expression of Ti oa d is of different forms in different 
applications. For example, when the engine is connected to 
a centrifugal compressor or a fan, the external load for the engine 
may result in a resisting torque proportional to the square of its 
rotational speed [25]. Therefore, in this study the load torque is 
assigned to be proportional to the square of the engine’s rotational 
speed with a load coefficient q oa d. 

In this study, the instantaneous gas forces exerted on the piston and 
the displacer are determined from the gas pressures in the expansion 
and the compression chambers with the help of the thermodynamic 
model [23]. On the other hand, the present dynamic model is used to 
calculate the instantaneous rotation speed of the engine (w) at any 
instant. The forces and torques exerted on the components and links of 
the engine at any instant are calculated by solving Eqs. (8)-(16) 
simultaneously as the gas pressures in the chambers are obtained 
from the thermodynamic model and the mass inertia, the friction 
force, and the external load are given. The instantaneous rotation 
speed of the engine w is determined by Eq. (18), and then the value is 
used to determine the new positions of the piston and the displacer at 
the consecutive time step, and then the gas pressures in the expansion 
and the compression chambers can be updated by the thermody¬ 
namic model again. Detailed information regarding the thermody¬ 
namic model has been provided in Refs. [23,24]; therefore, no further 
information will be provided herein to save space. 

Table 1 displays the geometrical and operating parameters of 
the baseline case. In the parametric study stated later, the variables 


Table 1 

Geometrical and operating variables of the baseline case. 


Variable 

Value 

Rd{ m) 

0.0036 

L(m) 

0.042 

h (m) 

0.018 

h (m) 

0.018 

h (m) 

0.018 

U (m) 

0.018 

n (m) 

0.020 

r 2 (m) 

0.0205 

Id (m) 

0.07946 

Ldt{ m) 

0.16374 

Lpt (m) 

0.05093 

Mm) 

0.158 

mi (kg) 

0.053 

m 2 (kg) 

0.028 

m 3 (kg) 

0.067 

m 4 (kg) 

0.067 

m 6 (kg) 

0.011 

P (kPa) 

101.3 

Th (I<) 

1000 

T c (I<) 

300 

Working gas 

air 

m a ir ( k §) 

3 x 10- 5 

£eff 

0.68 

Rn (°C/kW) 

19000. 

R t2 (°C/kW) 

5000. 

0 (degree) 

90.198 

Compression ratio 

1.7420 

Stroke (m) 

0.01 

Dead volume (m 3 ) 

2.93 x 10 6 

m 7 (kg) 

0.022 

m 8 (kg) 

0.022 

m 9 (kg) 

0.022 

m 10 (kg) 

0.022 


of the baseline case are changed to evaluate the effects of these 
parameters on the performance of the engine in order to find out an 
optimal combination of these parameters. All the studied cases are 
labeled in Table 2, and the geometrical and operating conditions 
associated with the studied cases are also given in this table. 

3. Results and discussion 

3.1. Baseline case analysis 

Fig. 4 shows the variation in the rotational speed of the engine 
during the start-up period for the baseline case. In the baseline 
case, the heat source and heat sink temperatures are 1000 I< and 
300 K, respectively, and the initial charged pressure in the engine is 
101.3 kPa. This figure shows the data of the instantaneous and the 
cycle-average rotational speeds. The engine is started from an 
initial rotation speed of 200 rpm, and eventually the cycle-average 
rotation speed reaches 915 rpm in the steady-state regime. The 
fluctuation of the instantaneous engine speed from the cycle- 
average value is not greater than 30 rpm in this case. The defini¬ 
tion of the cycle-average rotation speed is 

t+tp 

Wave(t) = T J wdt (20) 

t 

When the steady-state regime is reached, one has da) ave /dt = 0 
and then a combination of Eqs.(18) and (19) leads to 

'fload = "tpower (21a) 

or 
2 

w ave = T power/Ci oac | 


(21b) 
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Table 2 

Studied cases and their steady-state rotational speed, output power, and thermal efficiency at point of maximum power output. 


Case 

P(kPa) 

Th (K) 

Rd (m) 

G(m) 

Id (m) 

0 (degree) 

Compression ratio 

Stroke (m) 

Dead volume (c.c.) 

Values at point of maximum 
power output 

Wave (rpm) W (W) 7] 

1-1 

50. 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

912 

11.4 

0.4 

1-2 

101.3 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

789 

14.5 

0.435 

1-3 

200 . 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

650 

16.7 

0.45 

2-1 

101.3 

900. 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

742 

12.08 

0.43 

2-2 

101.3 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

789 

14.5 

0.435 

2-3 

101.3 

1100 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

826 

16.9 

0.44 

3-1 

101.3 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

789 

14.5 

0.435 

3-2 

101.3 

1000 . 

0.0036 

0.00005 

0.077 

90.198 

1.63 

0.01 

6.18 

748 

12.21 

0.36 

3-3 

101.3 

1000 . 

0.0036 

0.00005 

0.0754 

90.198 

1.57 

0.01 

8.29 

727 

11.12 

0.328 

4-1 

101.3 

1000 . 

0.002 

0.00005 

0.07946 

85.5 

1.35 

0.005 

6.39 

1008 

7.27 

0.21 

4-2 

101.3 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

789 

14.5 

0.435 

4-3 

101.3 

1000 . 

0.0043 

0.00005 

0.07946 

93.7 

1.96 

0.012 

1.14 

839 

20.7 

0.57 

5-1 

101.3 

1000 . 

0.0036 

0.00005 

0.07946 

90.198 

1.74 

0.01 

2.93 

789 

10.5 

0.435 

5-2 

101.3 

1000 . 

0.0036 

0.00004 

0.07946 

90.198 

1.79 

0.01 

2.93 

581 

12.43 

0.49 

5-3 

101.3 

1000 . 

0.0036 

0.00003 

0.07946 

90.198 

1.84 

0.01 

2.93 

467 

13.9 

0.51 


Since q oa d is a constant, it is expected that in the steady-state 
regime the cycle-average rotational speed is independent of the 
initial rotational speed. Fig. 5 reflects this expectation. Shown in 
this figure are the variations in the instantaneous and the cycle- 
average rotational speeds starting from different initial engine 
speeds, 150 rpm, 600 rpm, and 1200 rpm, for the baseline case. It is 
observed the initial rotational speed indeed has subtle influence on 
the variation in rotational speed in the start-up period; however, 
after the start-up period all different initial rotational speeds lead to 
the same final rotational speed in the steady-state regime. 

It is important to note that only an initial rotational speed within 
a certain range can successfully start the engine and make the 
engine reach the steady-state regime. It is practically reasonable 
that the engine is brought to a stop if an improper initial rotational 
speed outside the range is selected. 

The flywheel moment of inertia ( 73 ) is regarded as one of the 
influential parameters that dominating the rotation of engine 
during the start-up period. The effects of the magnitude of the 
flywheel moment of inertia for the baseline case are shown in Fig. 6. 
In this figure, the flywheel moment of inertias is assigned to be 
within the range between 4 x 10 -5 and 2 x 10 -4 kg m 2 . It is found 
that the fluctuation in the instantaneous rotation speed of the 
engine decreases with an increasing flywheel moment of inertia. 



This implies that a larger moment of inertia helps reduce the 
fluctuation in the instantaneous rotational speed. Meanwhile, in 
the cases with larger moment of inertia, longer time is needed for 
the engine to reach the steady-state regime. However, it is observed 
that the steady-state cycle-average rotational speed is not changed 
by altering the flywheel moment of inertia. This is attributed to the 
fact that in the steady-state regime doj ave /dt = 0 so that 
730:3 = l 3 (doj ave /dt) = 0. Therefore, as the steady-state regime is 
reached, the final cycle-average rotational speed is not affected and 
approaches to the same value as the moment inertia is varied 
within the range considered in the figure. 

In Fig. 7, the variation in cycle-average rotation speed at different 
load torques Ti oa d for the baseline case is displayed. The load torque is 
changed by adjusting the value of the load coefficient q oa d. In the 
cases shown in this figure, the values of ci oa d are set to be 9 x 10 -8 , 
3 x 1CT 7 , and 5 x 1CT 7 . In this study, the |3-type Stirling engine with 
rhombic-drive mechanism shown in Fig. 1, Model 10ST1, is used to 
measure the friction loss of the engine to evaluate q oa d. However, it is 
expected that the magnitude of cioad should be dependent on 
individual cases. Therefore, the test value is then extended to be 
9 x 1CT 8 , 3 x 1CT 7 , or 5 x 1CT 7 for a parametric study. 



Fig. 5. Instantaneous and cycle-average engine speeds vs. time, for baseline case 
starting from different initial speeds. 
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Fig. 6. Instantaneous and cycle-average engine speeds vs. time, for baseline case with 
different flywheel moment of inertias. 



The external load torque is applied from the very beginning, and 
thus the engine speed is decreased due to the existence of the load. 
Since the load torque is assumed to be proportional to the square of 
its rotational speed, as the engine rotation speed increases gradu¬ 
ally from the initial rotation speed to a higher speed, eventually the 
load torque is elevated to become equal to the power torque (T p0W er) 
in the steady-state regime. 

3.2. Parametric study 

Next, the effects of geometrical and operating parameters on the 
performance of the engine are evaluated. The performance of the 
engine is represented by the power output (W) and the thermal 
efficiency ( 77 ) of the engine. The present model is capable of pre¬ 
dicting the dependence of the performance of the engine on the 
rotation speed so that the performance curves can be plotted. The 
power output and the thermal efficiency are calculated after the 
steady-state regime is reached. The power output (W) is deter¬ 
mined by 



Fig. 7. Cycle-average engine speed at different load torques, for the baseline case. 


w — ^load ’ ^ave (22) 

and the thermal efficiency ( 77 ) of the engine is expressed by 


where Qi n is the input heat transfer rate, which can be determined 
by the thermodynamic model [23]. Fig. 8 shows the performance 
curves associated with the baseline case. The cycle-average engine 
speed ranges from 100 to 2500 rpm. It is interesting to find that the 
power output reaches its maximum at around 900 rpm, whereas 
the thermal efficiency reaches its maximum at around 420 rpm. 
That is, the rotation speed for maximum power output is different 
from that for maximum thermal efficiency. Therefore, it appears 
that for this particular case the optimal rotation speed is within 400 
and 1000 rpm, in which the magnitudes of both the power output 
and the thermal efficiency are relatively higher. 
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Fig. 9 conveys the effects of initial charged pressure on the 
performance of the engine, for Cases 1-1, 1-2, and 1-3 which are 
listed in Table 2. In this figure, the initial pressure is set to be 50, 
101.3, and 200 kPa, and the variation of the power output with 
the cycle-average rotation speed is plotted. It is observed that at 
Pjnitiai = 200 kPa (Case 1-3), the maximum power output is 
16.7 W at Wave = 650 rpm. The data for the steady-state rotational 
speed, the output power, and the thermal efficiency for these 
three cases at point of maximum power output are provided in 
Table 2. It is noted that the larger the initial pressure, the engine 
speed with maximum power output is lower, and the corre¬ 
sponding thermal efficiency and the maximum power output are 
higher. In this figure, it is observed that the power output is more 
sensitive to a change in initial pressure in low-speed regime than 
in high-speed regime. 

Fig. 10 shows the effects of the heat source temperature on the 
power output, for Cases 2-1, 2-2, and 2-3 in which the heat source 
temperatures are assigned to be 900,1000, and 1100 K, respectively. 
For the curve for Th = 1100 K, the maximum power output is 16.9 W 
at Wave = 826 rpm. The maximum power output is significantly 
increased as the heat source temperature is elevated. In addition, 
according to the data provided in Table 2, the higher the heat source 
temperature, the higher is the engine’s thermal efficiency. Contrary 
to Fig. 9, in this figure it is observed that the power output is more 
sensitive to a change in heat source temperature in high-speed 
regime than in low-speed regime. Therefore, an increase in the 
heat source temperature would be effective if the engine is 
designed to operate at higher speed. 

Effects of displacer length on power output are shown in Fig. 11. 
It is noted that as the stroke of the piston and the phase angle are 
fixed, a longer displacer means a higher compression ratio and 
a smaller dead volume. In this figure, the length of the displacer is 
set to be 0.0754, 0.077, and 0.07946 m, for Cases 3-1, 3-2, and 3-3. 
As shown in Table 2, the corresponding compression ratios are 1.57, 
1.63, and 1.74 and dead volumes 8.29,6.18, and 2.93 cc, respectively. 
It is expected that the case with a higher compression and a smaller 
dead volume is accompanied with a higher performance. In this 
figure, the results presented agree with this expectation. It is found 
that at Id = 0.07946 m, the maximum power output of 14.5 W is 




attained at 789 rpm. Also indicated in Table 2, the corresponding 
thermal efficiency is 0.435. The maximum power output and its 
corresponding thermal efficiency are reduced as the length of the 
displacer is decreased. For the case at Id = 0.0754 m, the maximum 
power output is only 11.12 W and the corresponding thermal effi¬ 
ciency is only 0.328. 

The effects of the offset distance from gear centers to joints 
between rhomboid and gears (Rd) could be very subtle since an 
increase in Rd results in increases in the compression ratio and the 
stroke of the piston as well as a decrease in the dead volume. 
Therefore, one may expect that an increase of Rd might generally 
improve the performance. However, since the phase angle will also 
be altered due to a change in Rd , the influence of this parameter 
becomes involved. Fig. 12 conveys the influence of Rd on the 
performance of the engine. In this figure, the cases at Rd = 0.002, 
0.0036, and 0.0043 m, for Cases 4-1, 4-2, and 4-3, respectively, are 



Fig. 12. Effects of Rd on power output, for Cases 4-1, 4-2, 4-3. 
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Fig. 13. Effects of regenerative gap size on power output, for Cases 5-1, 5-2, 5-3. 


investigated. It is found that in most cases, the power output 
increases with Rd. However, a changeover is found as the rotational 
speed is greater than a critical value located between 2100 and 
2200 rpm. In the higher rotational speed regime, larger is negative 
to the power output. That is, in this regime the larger the value of R& 
the lower is the power output. The reason for the existence of the 
changeover is still not clear and open to discussion at this moment. 
The change in the phase angle may play an important role on the 
changeover of the performance, and it is worth further investigation. 
However, with the rhombic-drive mechanism, the phase angle is not 
explicit and must be derived in terms of other geometrical param¬ 
eters of the mechanism. It is difficult to evaluate it sole effects by 
keeping other geometrical parameters fixed. 

As with Cases 5-1, 5-2, and 5-3 listed in Table 2, the varied 
displacer radius leads to different regenerative gap sizes (G) and 
compression ratios of the engine while the piston stroke, the 
clearance volume, and the phase angle are fixed. In Fig. 13, the 
power outputs as functions of cycle-average rotational speed 
associated with the three cases are plotted. The cases are at 
G = 0.00003, 0.00004, and 0.00005 m, which corresponding to the 
displacer radius of 0.00202,0.00201, and 0.002 m, respectively. The 
maximum power output is found to be 10.5 W at 789 rpm for 
G = 0.00005 m. When the gap size is decreased, the magnitudes of 
the maximum power output and the corresponding thermal effi¬ 
ciency of the engine are increased. In the case with G = 0.00003 m, 
the optimum power output of 13.9 W is reached at a lower speed of 
467 rpm. However, a decrease in the gap size is not favorable in the 
high-speed regime. Furthermore, a decrease in the gap size 
narrows the operating range of the rotational speed. For example, 
at G = 0.00005 m, the operating rotational speed ranges from 
approximately 100-2150 rpm, whereas at G = 0.0003 m, the 
engine can only operate between 100 and 1220 rpm. 

4. Concluding remarks 

A dynamic model for the beta-type Stirling engines with 
rhombic-drive mechanism has been developed and connected to 
the thermodynamic model [23] for dynamic simulation of the 
engine under different operating and geometrical conditions. 

Simulation of the baseline case shows that the instantaneous 
rotational speed of the engine always exhibits a fluctuating 


variation feature; however, the fluctuation of the instantaneous 
engine speed from the cycle-average value is not greater than 
30 rpm. The initial rotational speed obviously affects the variation 
in rotational speed in the start-up period; however, all different 
initial rotational speeds lead to the same final rotational speed in 
the steady-state regime. 

A larger moment of inertia helps reduce the fluctuation in the 
instantaneous rotational speed. Meanwhile, it is observed that the 
steady-state rotational speed is not changed by altering the 
flywheel moment of inertia. As the external load torque is applied, 
the engine speed is slowed down due to the existence of the load. 

Results of the parametric study show that as the initial pressure 
or the heat source temperature is increased, the engine speed with 
maximum power output is lowered and the corresponding thermal 
efficiency and the maximum power output are elevated. The power 
output is more sensitive to a change in initial pressure in the low- 
speed regime, but it is more sensitive to a change in heat source 
temperature in the high-speed regime. 

The maximum power output and its corresponding thermal 
efficiency increases with the length of the displacer or the offset 
distance from gear centers to joints between rhomboid and gears 
{Rd). However, in the higher rotational speed regime, larger Rd is 
negative to the power output, and thus in this regime a larger 
may reduce the power output. 

In addition, when the gap size is decreased, the magnitudes of 
the maximum power output and the corresponding thermal effi¬ 
ciency of the engine are increased; however, a decrease in the gap 
size narrows the operating range of the rotational speed. 
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